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Abstract. We present here a model for simulating the panchromatic spectral energy distribution of 
galaxies, which aims to be a complete tool to study the complex multi-wavelength picture of the 
universe. The model take into account all important components that concur to the SED of galaxies 
at wavelengths from X-rays to the radio. We review the modeling of each component and provide 
several applications, interpreting observations of galaxy of different types at all the wavelengths. 



INTRODUCTION 

The direct study of the process of formation and evolution of galaxies is now a realistic 
aim of observational cosmology. Thanks to facilities such as the HST, ISO, Chandra, 
XMM-Newton, SCUBA etc, now the theories of galaxy formation can be directly tested 
against multiwavelength observations of galaxies at different cosmological epochs. 

A panchromatic approach is critical both to obtain diagnostics of the various physical 
phenomena as well as to estimate the amount of reprocessing, due to the ISM, of 
the radiation produced by primary processes. In particular, a proper treatment of dust 
attenuation is compulsory when estimating the SFR from rest-frame UV luminosities 
(e.g. the cosmic SFR from Lyman Break Galaxies). This is particularly true if we 
consider that star forming regions are particularly dense and dusty environments. In 
a galaxy, dust reprocessing affects more the youngest stellar population, and the ensuing 
picture is much more complex than simple attenuation models frequently used. On the 
other hand, thermal emission from dust in these regions provides an alternative way to 
determine the star formation activity. 

To interpret the presently available observations, galaxy evolution models should 
include a careful treatment of the above mentioned effects. In this contribution, we 
review our population synthesis model GRASIL', which predicts time dependent SEDs 
of galaxies from X-rays to radio, including state-of-the-art treatment of dust and R^H 
reprocessing [1, 2], production of radio photons by thermal and non-thermal processes 
[3], nebular lines emission [4], molecular line emission, and X-rays expected from stellar 
populations. 



The code is available at http : / / adlibitum . oat. ts. astro, it /silva/def ault . html or 
http : / / web . pd .astro. it /granato/grasil/grasil. html 



DUST REPROCESSING 



One of the most important distinctive features of GRASIL is that it has included, for 
the first time, the effect of age selective extinction of stellar populations (younger stellar 
generations are more affected by dust obscuration), mainly (but not only, see below) due 
to the fact that stars form in a denser than average environment, the molecular clouds 
(MCs), and progressively get rid of them. 

Once the star formation and chemical enrichment histories of a galaxy are given, 
GRASIL computes the interaction between the stellar radiation and dust using a rela- 
tively realistic and flexible geometry for both stars and dust. In general, the system is 
described as a superposition of an exponential disk component and a bulge component, 
the latter modeled by a King profile. The ISM is divided into two phases, the diffuse 
ISM, corresponding to cirrus dust, and the much denser molecular clouds (MCs). Also, 
the fact that new stars are bom inside MCs and progressively get rid of them (either 
because they escape or because the clouds are destroyed) is taken into account. This 
process leads to an age selective extinction of stellar generations, in the sense that the 
younger the stars, the more they are affected by dust obscuration in MCs. This is de- 
scribed in GRASIL assuming that the fraction of starlight radiated inside the clouds by 
stars is a function of the star age. In practice, if ?esc is the timescale for the process, 
100% of the stars younger than tesc are considered to radiate inside the MCs, and this 
percentage goes linearly to 0% in 2fesc- The timescale ?esc is a fundamental parameter, 
which seems to be quite longer in starbursts than in normal disk-like galaxies [1]. 

The starting input for GRASIL is the history of star formation and chemical enrich- 
ment of the system. This is computed by some external code, which can result from a 
complex scenario for the formation of galaxies in a cosmological context (e.g. [5, 6, 7]), 
or a standard chemical evolution model for the formation of a single galaxy. Here, as in 
Silva et al. [1], we follow the latter approach. 

The total gas mass (diffuse+MCs) of the galaxy at time Tq is given by the chemical 
evolution model. The relatively fraction of molecular gas is a free parameter of the code, 
/mc- The total molecular mass, Mmc is then subdivided into spherical clouds of mass 
and radius, mmc and rmg. Then, the radiative transfer of starlight through the MCs and 
diffuse ISM. The original version of GRASIL presented in Silva et al. [1] used only a 
single population of MCs all with the same mass and radius, the code has been later 
updated to treat many populations, each one characterized by its mass mmc,!, radius r^c,i 
and total mass in the population Mmc,! (Silva, 1999, PhD thesis [8]). 

The ratio between gas mass and dust mass 5 is usually assumed to scale linearly 
with the metallicity of the residual gas. This quantity times m^c/r^^. determine the 
optical depth of the clouds. The predicted SED depend on mmc and rmc only through 
the combination mmc/^'mc which is the true free parameter. 

For the dust composition a mixture of graphite and silicate grains and R\H is adopted, 
with size distributions triggered to match the extinction and emissivity properties of the 
local ISM (for more details see Silva et al. [1]). However the code has the flexibility to 
modify the dust mixture, and to differentiate between MCs and cirrus. 

The SSPs included in GRASIL are based on the Padova stellar models, and cover a 
large range of ages and metallicities. Starlight reprocessing from dust in the envelopes 
of AGB stars is included directly into the SSPs, as described by Bressan et al. [9]. 



RADIO EMISSION 



The existence of a correlation between FIR and Radio emission is locally well estab- 
lished over a significant range of luminosity, from normal spirals to the most extreme 
ULIRGs. At 1.49 GHz, ^1.49 = log (Ffir/Fi. 49GHz /vfir) - 2.35±0.2 [10]. This correla- 
tion suggests that non-thermal (NT) emission is related to the recent star formation and 
the most likely mechanism is synchrotron emission from relativistic electrons acceler- 
ated into the shocked interstellar medium, following core collapsed supernova (CCSN) 
explosions [1 1]. Its small scatter testifies the universal proportions with which energy is 
radiated away at IR and radio wavelengths. Its validity has been recently confirmed up 
to redshift ~ 1.3 [12] and it is widely extrapolated much beyond, to estimate the redshift 
of more distant objects (e.g. [13]). 

In order to increase the diagnostic capabilities of GRASIL, Bressan, Silva & Granato 
[3] have extended its domain into the radio wavelengths. In brief, two main components 
were considered, the thermal emission from H II regions and synchrotron emission from 
relativistic electrons. The thermal component (Ly), is proportional to the number of 
hydrogen ionizing photons (Qa) derived from the adopted SSPs: 

(1) 

The non-thermal (Ly ^) component has been calibrated with the NT emission and the 
CCSN rate (rccsN) of our Galaxy, after accounting for the small contribution, ~6%, of 
radio supernova remnants (E^^^): 

Ly oc[t V +E V jrccSN (2) 

Though the calibration with our galaxy is quite uncertain, with these assumptions Bres- 
san et al. [3] and Panuzzo et al. [4] were able to reproduce well the FIR/Radio correlation 
of normal star forming galaxies, namely ^1.49 ~2.35. Bressan et al. [3] have also shown 
that, due to different spectral index of thermal and NT radio emission and the delay of 
CCSN explosions, radio emission can be used to analyze obscured starbursts with a time 
resolution of a few tens of Myr, unreachable with other star formation indicators. This 
is shown in Figure 1, from Prouton et al. [14]. A systematic study of compact ULIRGs 
from NIR to radio wavelengths is underway (Clemens et al. in preparation). This will 
provide important clues on the contribution of the AGN in ULIRGs and on the technique 
of photometric redshift determination of high-z obscured sources. 



MOLECULAR LINES 

The emission from molecular lines is used to diagnose the physical conditions within 

star forming molecular clouds. In fact, emission lines such as those from CS, HCN 
and HCO+ are excited only in high density gas, while others (e.g. ^^CO) trace more 
typical regions in the cold ISM. The study of the physical conditions of the star for- 
mation regions, the star formation efficiencies, the quantification of the contamination 
of molecular lines to the sub-mm continuum, and the possibility of obtaining spectro- 
scopic redshifts of distant galaxies using molecular lines are some of the challenges of 
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FIGURE 1. NIR to radio SED fits of two compact ULIRGs, IRAS08572+3915 and CGCG436-30. 
IRAS08572+3915 has a "flat" radio slope, due to free-free emission from HII regions. This is typical of 
"young" starburst. Its 1 .4GHz flux is attenuated by free-free absorption. The radio SED of CGCG436-30 
is already dominated by synchrotron emission from CCSNs, a feature typical of older starbursts. Upper 
solid line is the SED fit, dashed line is the contribution from the galaxy (starburst plus disk), dotted line is 
the contribution from the AGN and the lower solid line is the contribution of the disk alone. See Prouton 
et al. [14] for details. 



the next future millimeter astronomy. Therefore, the joint modeling of molecular and 
dust emissions is essential to investigate the global relationships between gas, dust and 
star formation. 

We have developed a method that combines GRASIL and a molecular lines emission 
model, based on LVG approximation (Vega et al. in preparation). Our LVG code (see 
de Jong et al. [15] for the mathematical formulation) takes into account coUisional 
excitation by H2 molecules, radiative trapping and radiative excitation by the cosmic 
microwave background. The adjustable LVG parameters are: the numerical density of 
colliders, nn^, the gas kinetic temperature, Tk, and the parameter A = X^oi'^mc/Ay in pc 
km^^, where X^^oi is the relative abundance of a given molecule to H2. Our LVG code at 
present treats ^^CO, ^^CO, C^^O, CS, HCN, HNC and HCO+. 

The main advantage of combining the LVG code with GRASIL for the analysis of 
molecular emission in galaxies is that, by fitting the galaxy SED from the far UV to the 
radio bands, GRASIL can provide realistic estimates of the SFR, fraction of molecular 
gas, gas to dust ratio, dust temperature, and the average molecular cloud structure. 

The shape of the predicted IR SED depends mainly on the optical depth of the 
clouds, Tmc °^ S ■"^mc/'"mc- After is Set from the fit to the continuum, we obtain 
the molecular gas density = O.Vmmc/ (A-nr^^-mnj^/S) corresponding to different mmc 
(within reasonable values), assuming that most of the molecular emission comes from 
the same clouds with ongoing star formation. We apply our LVG model to the clouds 
for each value of mmc (or nn^), adjusting Tk and A to fit the observed intensity of 



40 



50 



20 



10 



NGC7469 



R1056+24HII 
HII + Liner 



NGC57INS 05052' 
:NGC61S1,-S|'' ffll 
HII 



....................... 

10.0 10.5 11.0 11.5 12.0 12.5 
Log (Lfir) 

FIGURE 2. Kinetic temperature versus FIR luminosity. The line represents the linear fit for the starburst 
galaxies. 



the molecular lines lines. The best fit model gives the best estimate of the excitation 
conditions of molecular gas. 

As an example, we show some results of applying our GRASIL+LVG model to six 
galaxies, four LIRGs and two ULIRGs, for which both continuum and molecular line 
data are available. From our preliminary analysis, there are hints that Tk is tightly 
correlated to the FIR (and Radio) emission over a wide range of luminosity (Figure 
2). The only galaxy out of this correlations is a Syl.2. A likely explanation is that, 
in starburst galaxies, the molecular excitation is driven by processes related to star 
formation. In the case of the AGN, the molecular gas seems to be "overheated", probably 
by the central engine. The mean molecular densities found for these galaxies are < 10^'^ 
cm^^ and the bulk of CO emission is coming from a low density medium. However, 
observations of HCN suggest the presence of a much denser medium, thus we are 
working in a two-zone model: dense clouds embedded in a low density medium. 



X-RAY EMISSION FROM GALAXIES 

The X-ray range has recently become very effective to study galaxy evolution, thanks 
to the high resolution, low flux levels and large spectral range reached by Chandra and 
XMM-Newton, that allow to probe the faint X-ray sources of star forming galaxies. In 
several works a strong correlation between the X-ray and the FIR or radio luminosity 
has been reported for star forming galaxies (e.g. [16, 17, 18, 19]). Grimm et al. [20], by 
comparing the luminosity functions of X-ray point sources (XLF) of several galaxies, 
conclude that the shape of the XLF being almost universal. Therefore, the X-ray emis- 
sion can be considered another independent indicator of star formation rate and history, 
to exploit together with the radio, IR and UV indicators, in order to set new constraints 
for a multi- wavelength tackling of galaxy evolution studies. To complete our self con- 
sistent modeling of the SED of galaxies, we have extended GRASIL to the X-ray range 
([21] and Silva et al. in preparation). 
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FIGURE 3. X-ray SED for M82, obtained from the SF history of the model fitting the UV to radio 
data. In the first plot we have used the X-ray predicted by including the SSPs of all the metallicities (from 
Z — 0.0004 to 0.05), while in second plot we have used only the Z — 0.008 SSP. In the latter case, even 
though the wind mass loss rate is lower, the higher BH remnant masses determine much higher accretion 
luminosities. Note however that the stellar inputs (initial-final mass relation, its dependence on metallicity, 
initial mass for the formation of a BH...) are very uncertain. Here an absorbing column density Nh = 10^^ 
cm^^ has been applied to the model SED. The observed SED (dashed line from 2 to 10 kev) is from 
Griffiths et al. [25]. 



We consider the following sources (see [22, 23]): (1) Binaries in which the donor is 
a core H-buming OB star that fills its Roche lobe, and the primary is a black hole (BH) 
or a Neutron Star (NS) (hereafter RBH and RNS); (2) Binaries in which the donor is a 
shell H-burning OB star undergoing a phase of strong stellar wind, and the primary is 
a BH or a NS (hereafter WBH and WNS); (3) Pulsars (PLS); (4) Supernova remnants 
(SNRs); (5) low-mass X-rays binaries (LMXRBs). 

The X-ray luminosity from the different components is directly included in the library 
of Simple Stellar Population. We proceed in the following way: (i) at each age of a 
stellar population we know the number of stars in the different evolutionary phases or 
dead remnants, therefore we compute how many stars of main sequence or post main 
sequence combine to form binary systems with BH or NS as a primary; (ii) then each 
source is assigned the appropriate Lx, therefore we obtain the total Lx of a population 
of objects of each kind; (iii) we distribute this energy according to suitable (observed 
average) spectral shapes. 

The inputs needed for the computation are the following (see [24, 22]): the IMF for 
single and primary stars, the distribution of the initial mass ratio M2/M1, the binary 
frequency at birth (> 0.5), the distribution over the periods. Of course the results depend 
also on assumptions on stellar evolution (mainly mass loss rate, remnant masses, initial 
masses for BH formation, dependence on metallicity). As for LMXRB, they are linked 
to the old stellar populations, i.e. to the mass content of galaxies [20]. We adopt the ratio 
Lx/M^ evaluated for LMXRBs in our Galaxy by Grimm et al. (5 x 10^^ erg/s/M©) and 
assume this ratio is the same in all galaxies. 



Once we have the SED of the populations of different X-ray sources as a function of 
age and metallicity, the integrated X-ray SED for a model galaxy of age Tq is given by: 



where \if{t) and Z{t) are the star formation rate and metallicity enrichment history for the 
model, i refers to the different X-ray sources, and fj{£) is the spectral shape associated 
to source of kind j. In Fig. 3 we show the X-ray SED resulting from our model for 
M82. The star formation and metallicity evolution for M82 are obtained by fitting its 
UV to radio SED [1, 3]. The plots show the important effect on assumptions on stellar 
evolution and metallicity. Work is in progress to explore these effects, also connected 
with the presence of the so called Ultra-Luminous X-ray sources (ULX) in star forming 
galaxies (e.g. [26]). These sources are interpreted as intermediate mass BHs (M > 100 
M0) or to specific effects of the accretion onto normal BHs [27, 28] with a violation of 
the strict Eddington limit. 



Emission lines from Hll regions are a powerful diagnostic of the instantaneous SFR 
because the luminosity of H recombination lines is proportional to the number of mas- 
sive, short-living stars. Moreover, emission lines are widely used to estimate the dust 
attenuation, e.g. with the Balmer Decrement. Nebular lines are also used to estimate the 
metallicity of the gas and, finally, to infer the nature and age of the ionizing source. For 
all the above reasons, we implemented the computation of nebular lines in GRASIL (see 
Panuzzo et al. [4] for details). 

To compute the line emission intensities in a population synthesis model, one should 
consider the spectrum of the ionizing source, provided by the recent star formation 
history, and use a photoionization code with suitable values of the gas parameters. 
It becomes particularly time consuming in applications requiring a large number of 
models. Our approach has been to pick out the physical parameters which actually affect 
the emission properties of H II regions. 

The emission spectrum from a single H II region depends on two main ingredients: 
the SED of the ionizing star cluster (determined by the IMF and the total mass, age 
and metallicity of the cluster, as well as by the adopted model atmospheres) and the 
properties of the excited gas (the density, the chemical composition and the geometry). 

We found that the emission line spectrum of an H ll region with fixed gas properties 
is described with reasonable precision by only three quantities: the number of ionizing 
photons per second for Hi, He I and O ii ((2h> Qho and Qq) emitted by the ionizing 
source. We used as ionizing SEDs analytical spectra composed by piece-wise black- 
bodies, similar to young SSP spectra, with temperatures given as functions of Q^, Qyi^, 
and Qq. We shown that our analytical spectra produce the same emission line spectra 
of SSP with the same values of (2h> !2He> and Qq, within a reasonable accuracy. Thus, 
we computed a library of photoionization models as a function of Qy{, Qne^ and Qq. 
When estimating the actual line emission due to a given stellar population, the spec- 
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FIGURE 4. The A^y derived from the ratio Ffir/^uv vs. Ahh from Balmer decrement. Filled circles 
show Buat et al. [29] sample, compared with models (different symbols connect by lines). The Auv = 
1.6AHa is the observed correlation for starbursts. See Panuzzo et al. [4] for details. 

trophotometric code computes Qn, Qn^, and Qo from the corresponding SED, and then 
interpolates the value from the above library. 

The library was computed for a grid of values of Qn, Qhs/Qh, and Qo/Que, for 
different assumptions on the gas density, metallicity and filling factor £. This method 
allows to get rid of the particular SSP model and IMF. Moreover it can be implemented 
in every population synthesis code. The library of photoionization models is available to 
the astronomical community from a web site^. 

A simple application of our population synthesis with nebular emission and dust 
processing was to study different methods to estimate the attenuation in normal spirals 
[4]. UV attenuation (Auv) can be estimated by i^piR/^uv ratio (see [30, 31]) while the 
attenuation suffered by Ha (Aua) can be derived from the Balmer decrement. Observing 
a sample of normal spiral galaxies, Buat et al. [29] found that Auv and Ahk do not 
correlate, contrary to what found in starbursts galaxies [32, 29]. 

Our models [4] (see Fig. 4) reproduce location of the observed galaxies and confirm 
a real lack of correlation between the Auv from Ffir/Fuv and Ana. The dispersion is 
a clear effect of the age selective extinction. In fact massive ionizing stars have a short 
lifetime so they can spend all their life inside dusty environments (i.e. MCs) while UV 
emitting stars that have a longer lifetime can radiate from a less dusty medium. 



PAH MODELING 

The mid infrared (MIR) spectra of many dusty galactic and extragalactic objects show 
strong emission features most commonly ascribed to aromatic C-C and C-H vibrations 
in large planar Polycyclic Aromatic Hydrocarbons (PAH) molecules (e.g. [33]). 

Through a comprehensive analysis of the MIR to FIR spectra of Galactic star forming 



http : / / web . pd .astro . it /panuzzo /hi i 



regions, normal and starburst galaxies and AGN, Peeters et al. [34] have shown that the 
SED of galaxies is very similar to that of Galactic exposed PDRs, and that PAHs are a 
better tracer of B stars rather than by massive star formation (O stars). This has been 
found also by Tacconi-Garman et al. [35] by direct high resolution 3jUm imaging of two 
nearby starbursts (NGC253 and NGC1808): they find no spacial coincidence between 
the detailed distribution of PAH emission and the locations of most recent star formation. 

On the theoretical side, the new observations have allowed to better define the model- 
ing of PAHs (e.g. [36, 37]). In order to interpret at best the spectra of galaxies, a detailed 
modeling that includes the complex interaction between stars and dust is needed. To this 
aim we have updated GRASIL by including a better treatment of PAH emission bands. 
The previous version of GRASIL [1,8] included PAH modeling based on pre-ISO data: 
the five 3 to 11 /im PAH bands were computed following mainly Xu & De Zotti [38]. 
We have now updated our treatment according to the PAH model by Li & Draine [36] 
(LDOl). In particular: (i) ISO observations have increased the number of known bands, 
therefore now in addition to the 3.3, 6.2, 7.7, 8.6, 1 1.3/im bands, we compute the 11.9, 
12.7, 16.4, 18.3, 21.2, and 23.1jUm bands; (ii) we adopt a Drude profiles for the fea- 
tures; (iii) the FWHM and emission cross section for the new bands are from LDOl. We 
have compared our updated PAH model with ISO spectra of galaxies by Lu et al. [39] in 
order to test weather the LDOl model, suited for the diffuse ISM, can be safely adopted 
in general for star forming galaxies. We searched continuum NIR to radio data for the 
galaxies of the Lu et al. sample, and fitted the continuum SED. We then made a detailed 
comparison of the corresponding model PAH with the Lu et al. MIR spectra. In all cases 
the predicted PAHs compare very well with the observation, the greatest disagreements 
are always within a factor of 1.5, but typically much less. 



CONCLUSIONS 

GRASIL is a realistic model for simulating detailed and panchromatic SEDs of galaxies. 
We are constantly working to improve GRASIL modeling in particular about PAH, X- 
rays, molecular line emission, and dust geometry. The model was successfully used in 
many works and in several different fields, also by researchers outside the developing 
team. GRASIL is available for the scientific community: everybody is welcome to 
download the code from our web site or to ask for the source. We are now working 
on a web interface: any researcher will be able to run GRASIL on our computers, and to 
get the result by e-mail. 
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